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Abstract 

In this paper, we study the role of boundary conditions on the optimal shape of a dyadic tree 
in which flows a Newtonian fluid. Our optimization problem consists in finding the shape of the 
tree that minimizes the viscous energy dissipated by the fluid with a constrained volume, under 
the assumption that the total flow of the fluid is conserved throughout the structure. These 
hypotheses model situations where a fluid is transported from a source towards a 3D domain 
into which the transport network also spans. Such situations could be encountered in organs 
like for instance the lungs and the vascular networks. 

Two fluid regimes are studied: (i) low flow regime (Poiseuille) in trees with an arbitrary 
number of generations using a matricial approach and (ii) non linear flow regime (Navier- Stokes, 
moderate regime with a Reynolds number 100) in trees of two generations using shape derivatives 
in an augmented Lagrangian algorithm coupled with a 2D/3D finite elements code to solve 
Navier-Stokes equations. It relies on the study of a finite dimensional optimization problem in 
the case (i) and on a standard shape optimization problem in the case (ii). We show that the 
behaviours of both regimes are very similar and that the optimal shape is highly dependent on 
the boundary conditions of the fluid applied at the leaves of the tree. 

Keywords: Shape optimization, Dyadic trees, Poiseuille's law, Fluid Mechanics, Stokes and 
Navier-Stokes systems. 
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Introduction and motivations 



Tree structures are very common means to transport a product between two regions of different 
scales. A fluid acting as a transporter often flows in such structures. The circulation of the 
fluid in such geometries can dissipate a lot of energy by viscous effects and the question of the 
optimization of the tree geometry arises. Important applications of this problem exist, like in 
industry, in the study of river basins, in water treatment, in medicine. Such structures are also 
often encountered in Biology and are the result of evolution. One can see roughly natural selection as 
an optimization process which adapts the organisms to their environment. Thus, a major question 
arises for biological systems: what are they optimized for? We know that an effect of natural 
selection is to minimize some complex cost functions relatively to a wide range of parameters. 
Although most of the time unknown, these parameters have probably various influences in term of 
amplitude, some stronger than others. Thus, if we are able, typically through a modeling work, 
to give hypotheses on what are the most influential parameters, to isolate them in a model and to 
determine their role if they were alone, then a simple comparison between the results of the model 
and the real biological system could indicate whether their role was truly important or not. In the 
case of biological networks such as lungs or vascular network, two cost parameters arise naturally: 
the viscous dissipated energy of the fluid in the network and the volume of the network. Indeed, 
these organs have to deal with energy dissipation due to air or blood circulation [33] and since 
they span in the geometry they have to feed, they cannot use too much volume. These organs can 
be subject to dysfunctions which are often consequences of an increase of their hydro dynamical 
resistance and thus of a loss of efficiency of the geometrical structures as transport systems. Hence, 
dysfunctions like asthma or heart attacks are often linked to an increase of the hydrodynamic 
resistance (energy cost spent for the circulation of the fluid) of the structure. Thus, researching the 
optimal shapes of tree structures could bring useful information. Previous studies have been made 
on this topic in the past like in [161 121 ESI [2Ql E] each on particular situations and applications. 

In this frame, the goal of this work is to determine the shapes of dyadic trees that would minimize 
the viscous energy of a Newtonian fluid under a volume constraint on the tree. As written upwards, 
such structures are good candidates for the modeling of mammals bronchial trees [20] and we will 
focus most particularly on this application. We study this problem for two regimes of flow. We 
begin with low flow regime (Poiseuille flow) using a matricial formulation of the problem as in 
[22\ 12 lj . We also study the non linear Navier-Stokes flow at a moderate Reynolds number (~ 100) 
which is however sufficient to exhibit inertial effects around the bifurcation. In this second case, 
we use a numerical shape minimization method based on shape derivatives [13]. We show that for 
both regime the optimal structure depends on the fluid boundary conditions that are applied at tree 
root and leaves. Indeed, under the assumption that the total viscous flow is constant throughout 
the tree, the optimal shape is very different according to the boundary conditions imposed at the 
leaves: Dirichlet conditions or strictly identical Neumann conditions lead to an optimal tree with 
particular relationships between the flow in a branch and its diameter; non identical Neumann 
conditions lead to a degenerated tree reduced to a tube and only one leaf remains accessible to the 
fluid from the root. Moreover the numerical simulations in the non linear case give precisely the 
geometry of the bifurcation. We focus here on the 3D case, however it is easy to see, with very few 
changes in the reasoning, that our results would also hold for the 2D case. 
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1 Terminology and notations 



In the following, we will call inlet of a dyadic tree either the open surface of the root of the tree or 
the root of the tree itself, depending on the context. Similarly, we will call outlets of a tree either 
the open surfaces of its leaves or its leaves themselves, also depending on the context. Mainly, we 
will refer to the open surfaces if we speak of boundary conditions and to the branches in the other 
cases. Note that this terminology does not mean necessarily that fluid is going in the tree through 
the inlet and out of the tree through the outlet. 

We will use the following notations throughout this section: 



[x], with x£l 

B T , where B is a matrix 

(ei, . . . , e n ), with n G N* 

(x, y), where x and y are two vectors with same 
length 

diagu, where u = (ui) ie ^ n j G R nxl 



the integer part of x 
the transpose of B 
the canonical basis of R 
the euclidean inner product 

the euclidian norm, induced by the inner prod- 
uct (., .} 

the diagonal matrix D = (dij)l<i,j<n such that 
di t i = Ui for all i G 



In the following, a family of real numbers (v £ )s€R* + is assimilated to a sequence since for each 
n G N*, one can choose e = —. 



Moreover, let us define the direct sum of two matrices. 

Definition 1. Let m and n be two nonzero integers. Let A\ G R™ x?1 and A2 G R mxm be two 
matrices. The direct sum of A\ and A2 is the matrix M in j^™+ m >™+ m defined by blocks by 



M = Ai e A 



A, 








A 2 



2 Models 

We introduce here all the models used in this 
section have been regrouped in Appendix lAl 



article. For the sake of clarity, all proofs of this 



2.1 Poiseuille's law 

We consider a viscous incompressible fluid whose dynamic viscosity is \i > 0. It flows in a steady 
and laminar state through a cylindrical rigid pipe whose length is L > and radius is R > 0. We 
impose a no-slip condition on the lateral boundary which means that the fluid "sticks" to the wall. 
We refer the inlet to and the outlet to 1. Pressures (po > and p\ > 0) at its openings are 
supposed to be uniform all over the section. The volumetric flow rate is chosen positive if the 
flow goes from section to section 1. 

The fluid behavior is ruled by the Navier-Stokes equations in the cylinder. The solution of these 
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equations with the previous conditions is characterized with velocity profiles that are parabolic on 
each section and a pressure that is constant on each section and decrease linearly along the axis 
of the pipe. This type of flow is known to verify Poiseuille's law. This law boils down to a linear 
relationship between the pressure drop po — p\ and the volumetric flow rate through the pipe, 
that is 

Po~Pi = r§ 

where 

8uL 



Therefore, if the pressure drop po — pi is positive, the flow goes from section to section 1. 

There is an exact correspondence between Poiseuille's law and Ohm's law if we match pressure 
drop with potential difference and volumetric flow rate with electric current. Thus by analogy, the 
proportionality constant r is called a hydrodynamic resistance. Finally, since the flow is incom- 
pressible, the flow rate is conserved through the pipe. 



2.2 Dyadic trees 

In this section, we will define the different mathematical tools needed in the sequel to manipulate 
tree structures in which flows a Poiseuille's fluid. In particular, we give a matricial relationship 
between the flow at the inlet of a dyadic tree and the pressures at the outlets. 

We consider the flow of an incompressible and viscous fluid whose dynamic viscosity is \i through 
a finite dyadic tree of N + 1 generations (N € N, N > 1). We recall that a new generation of this 
tree occurs when a bifurcation is created. Hence the root branch corresponds to generation 1 and 
the branches at the leaves of the tree corresponds to generation N + 1. Consequently, the tree has 
2^ outlets and 2 N+1 — 1 branches. 

We also define the level as a number associated to a generation, equal to 1 for the second 
generation and increased by 1 at each generation. Therefore, a tree with N + 1 generations has N 
levels. 

We do the distinction between the generation and the level in the tree, because it makes the 
indexation of the different variables of the tree easier. The levels will be denoted in all this section 
by the index letter i. 

We assume that our tree is composed of connected rigid cylindrical pipes in which the fluid 
obeys to Poiseuille's law. We call $ > the volumetric flow rate that enters the tree. Since the 
flow is incompressible, the flow rate is conserved through the pipes and at bifurcations. 

The sets Bn,i of couples of indexes that locate each branch of a given level % (i £ [l,iV]) is 

B N ,i = {(ij)\j€ii,n}. (i) 

Hence, i represents the level and j the position of the branch at this level. Thus, the set Bn of all 
indexes locating the pipes for the overall tree is 

Bn= (J B N , t (2) 
ie[i,7V] 
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This set does not include the root branch of the tree. This branch has a radius Rq > 0, a length 
Lq > 0, the pressure at its inlet is po > (pressure at the inlet of the tree) and is p\ > at its 
outlet. Thus, its hydrodynamic resistance is ro = According to Poiseuille's law, one has 

TTU 



Po-Pi = r $. 

For a pipe whose location in the tree is given by the couple € B^, we denote Rij > its 
radius and Li ,• > its length. Therefore, the hydrodynamic resistance of this pipe is rj = p^-. 

We use qij and pij > to define respectively, the volumetric flow rate through this pipe and the 
pressure at its outlet. The flow rate qij is chosen positive if the fluid in a pipe flows towards the 
pipes with a higher generation index. 

We consider now a pipe whose index is (i, j) E B/v.i- Since the tree is dyadic, the 

ie[i,iv-i] 

hydrodynamic resistances of its (two) daughter pipes are rj+i^j-i and r^i^j, which belong to the 
(i + l)-th level. Then, we define the reduction ratios Xj+i^j and Xj+i^j+i, that represents the 
change in the geometry of the pipes between the levels i and i + 1, by 

r i,j Lij ( Ri+l,2j-l \ 4 , Tij Lij ( Ri + \2j^ 

x i+h2 j-i = = ^ — and x i+h2 j = — = j — „ 

r i+l,2j-l ^i+l,2j-l \ J r «+l,2i L i+l,2j \ *H,j 

(3) 

Moreover, we assume an identical reduction ratio for the radius and the length between a mother 
branch and its daughter, thus 

( Ri+l,2j-l\ 3 , fRi+l,2j\ 3 (A , 

Xi+i,2j-i = ^ and x i+h2j = — ^ . (4) 



Figure Q] shows an example of a dyadic tree and of our notations in the special case where N is 
equal to 2. 




% = 1 
i = 2 



P2,4 



Figure 1: An example of dyadic tree with three generations (N = 2). Note that the flows q% can 
be either positive or negative £ Bn)- 
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Our goal is to establish a relationship between the 2 N pressures and the 2^ volumetric flow 
rates at the outlets of the tree. To go further, we need to be able to follow the fluid through paths 
in the tree. Therefore, we define the notions of path and subpath in the tree. 

Definition 2. Notion of path and subpath. 

1. A path H _^uj\ (with € Bn) is the set of couples of indexes of the i branches needed to 
link the root branch denoted by and the branch located by in the tree. It includes the 
branch referred to by (i,j) but not the root branch. More precisely, 

n 0-Ki,i) = {(!> m i) ,•••>(*- 1) m 2) , (i, mi)} , (5) 
where (mk)i<k<i denotes the sequence of positive integers defined by 

m 1 = j 



m k+1 



m k + l 



, Vfce [M-l]. 



2. Let € Bn and s G [0, ij. For a given path Hq^u^, the subpath Hq^u^s) is the set of 
couples of indexes of the m branches needed to link the root branch and a branch located at 
the s-th level following a part of the path n _j.(ij). More precisely, the subpath IIo_>(j j) (s) is 
the subset of n _ 5 .(j i j) defined by 

rr r a \ - / {(l,mi),...,(s,mi_ s+ i)} if s > 1 ( 
n o^)(«)-{ ifs = _ (6) 



We use this definition to do a change of variable. It will be very useful in the second part of 
this section devoted to the optimization of a given criterion with respect to the geometry of the 
tree represented by the 2 N+1 — 2 variables { x i,j}(i,j)eBN' From now on, we replace the variables 
Xij by the new ones £jj defined by 

V(i,j) G B N , £ it j = Yi x k,i- ( 7 ) 

(fc,i)en _>(i, 3 -) 

We will impose another constraint on the geometry: lengths and radii of the tree are assumed 
to decrease as we go along its levels. More precisely, 

V(iJ) € |J B Nji , max(^ + i )2 j-i,Ci+i,2i) < (8) 
»e[Mr-i] 

It has to be noticed that the map {%i,j}(i,j)eB N 1 — y i^i,j}(i.j)eB N 1S obviously a C 1 -diffeomorphism 
on the set of strictly positive real numbers so that it defines a change of variable. 
Moreover, represents the hydrodynamic resistance r^j of the pipe denoted by 
For instance, in the case N = 2 (see figure [rj, the path linking the inlet to the first outlet is 
IIq_>.(2,i) = {(1, 1); (2, 1)} and the geometric variables of theses branches are 

ro ( #i,i \ 3 c , r f R 2 , A 3 t 
^1,1 = = ' «,1 = and x 2,l = = -pr~ ) «,l = £1,1^1,2- 

n,i V R o J r 2,i V Ro J 
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Since the tree is dyadic, it is usual to compute the total volume V of the tree by summing the 
volume of each cylindrical branch, i.e. 

Volume = ttRILq + ^ ^R 2 i,j L i,j = vR^Lo f 1 + 

(ij)GSjv \ (i,j)eB N 




Notice that this volume is an approximation of that of the real tree, since it takes into account 
only the volumes of the cylinders composing the tree and not the volumes of the bifurcations. 

Now, let us define 

• the vector p containing all the pressures at the outlet of the tree, i.e. 

P=(.PNj)] ell , 2Nj e(K) 2Nxl , 

• the vector q containing all the volumetric flow rates at the leaves of the tree, i.e. 

Q = (<?Ar,j)J e [ 1)2 iV] € M 2 xl , 

• the vector £ representing the resistances of the tree, i.e. 

€ = (6,1,6,2, • • • ,6v,i, • • • ,6v, 2 ~) T e (K) (2N+1 ~ 2)xl - 

Definition 3. Given two positive integers i and j and their binary expressions 

CO CO 

i = a k2 k , j = ^ with (a k ,/3 k ) e {0, l} 2 , Vk e N, 

k=0 k=0 

we define v^j as 

Vij = min{/c € N | a\ = A, V7 > k}. 

Thanks to Poiseuille's law, we are able to establish a linear relationship between the vectors p 
and q. 

Proposition 1. Let N € N*. One has 

p u N -p = A N (£)q, (10) 

where 
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the symmetric matrix An{£) € R 2JVx2JV is called resistance matrix of the tree and is defined 
by 



N 



An(€) - ( a i,j)l<i,j<2*, ">ij 



(fc,OGn ^ {JV>i) (Af-j/ l _i, J _i) 



(11) 







mat G 



P 2 JV xl 



denotes the ones vector i.e. ujv = (!)•••) 1) T - 



Remark 1. n _ > (7v,i)(-^ — ^i-i.j-i) * s ^ e subpath corresponding to the intersection of the two paths 



This proposition is proved in appendix IA.1I A close result has already been proved in |22j . 
The quantity 

£(q,0 = q T A N (£)q (12) 

corresponds to the total viscous energy dissipated by the fluid in the tree during one second, see 
[221 EE], i.e. 



2™ „2 



2 JV ~ 1 \ 2 /^2 JV \ 2 

51,1 6,2 i^l^' j^l 



(13) 



Using the fact that the flow rate is conserved all along the pipe, we rewrite the dissipated energy 



ro 
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As a consequence, An(£) is a symmetric positive definite matrix. It justifies the following 
proposition. 

Proposition 2. The resistance matrix An(£) € M 2JVx2JV is invertible. 

Example 1. So as to have an idea of the structure of the matrix An($,), let us use an example 
of a tree with N = 2 levels (thus, with N + 1 = 3 generations, 2 N ~ 1 = 4 outlets and 2 N — 1 = 7 
branches). Its resistance matrix ^2(£) G M 4x4 is defined as follows: 



1 + — 1 + — + — 

^ 6,i ^ 6,i ^ 6,2 



V 



1 
1 

1 + 6,2 + 6,3 



1 + 



6,2 



1 1 1 + 7^ 1 + ^ + ?^ / 

6,2 6,3 6,4 / 

The tree corresponding to this example is drawn on the figure [0 



S 



2.3 Boundary conditions 



For biological systems like the vascular network or the lung, it is very difficult to know what are 
the correct boundary conditions since there exists complex feedback loops that are able to modify 
the forces applied by muscles relatively to the behaviour of physiological values. Typically, lungs 
ventilation is controlled not only by carbon dioxide and oxygen concentration in blood but also 
by blood pH and by mechanical sensors distributed along the tree [28]. Concerning boundary 
conditions in the frame of fluid mechanics, one can refer to [5] and [IT] and in the context of blood 
flow modeling to [27] . 

We will investigate the two major types of boundary conditions at the openings of the tree. 
We will either impose a flow (or parabolic velocity profile) which corresponds to a quantity of fluid 
that would circulate in the branch, or impose a pressure which corresponds to a force that would 
be applied to the opening of the branch. If flow rates are imposed at the outlets of the tree, then 
the pressure will be imposed at the inlet and vice versa so that the problem is well posed. Indeed, 
the two cases studied in this work are 

• 1st case: pressure is imposed at the inlet and flow rates at the outlets. 

• 2nd case: flow rate is imposed at the inlet and pressures at the outlets. 



In the 2nd case, the flows at outlets are not directly accessible. Thus, the following proposition 
gives the existence of q as the solution of a linear system. 

Proposition 3. Assume that the pressures at the outlets of the tree are fixed and the flow in the 
tree root equal to <I>. Then, the vector q G R 2 xl is the unique solution of the linear system 

M N (£)q = b N (14) 

where 

( {A N {i) Vl y \ ( -<p,vi) \ 



(A N (£)v 2N _ 1 ) T 



\ 



UN 



GM 2 x2 , b N 



/ 



D 2«xl 



(15) 



and for all i G [1, 2^ - lj, Vi = (0, . . . ,0, 1, —1,0, . . . ,0) T G IR 2iVxl , with 1 at the i-th position and 
— 1 at the (i + l)-th position. 

See appendix IA. 21 for the proof of this proposition. 



3 The optimization problems 

In this section, we will study the two finite-dimensional constrained optimization problems corre- 
sponding to two types of boundary conditions: 

• 1st case: pressure is imposed at the inlet and flows at the outlets, 



9 



• 2nd case: flow is imposed at the inlet and pressures at the outlets. 

We consider the same rigid dyadic tree as in section [2] with N +1 generations (N G N, N > 1) 
through which flows the fluid introduced previously. The tree is characterized by its geometry £ 
(£ G (M+) (2JV+1_2 ) xl ), its resistance matrix An(£) S M 2iV><2JV ) and the volumetric flow rates 

q (q G M 2 xl ) at its outlet. We recall that, at outlets, pressures are linked to flow rates with the 
relation p = A]y(£)q. 

We want to minimize the total viscous dissipated energy £ defined by (fT2j) with respect to the 
pair (<jf,£) under the constraints 

• that q is the flow vector at outlets of the tree, when the Poiseuille's model is considered. This 
condition can always be treated as a constraint, it takes however different forms depending 
on the case considered. 

• and that the total volume V of the tree is constrained, which from Q, writes 

1 + E ^ = A ( 16 ) 

with A = X T > 1. 

To make the admissible set of pairs (q, £) clear, we define the intermediate set g/\ as 

a/ K = || g (M;) (2iY+1 " 2)xl | £ verifies the conditions © and ([IB]) } . (17) 

3.1 1st case: pressures at the inlet and flow rates at the outlets are known 

The first case is much simpler than the second one (presented in section I3.2|) . The pressure is 
imposed at the inlet along with the flows vector q at the outlets, thus the optimization problem is 

where q G IR 2 ^ xl is assumed fixed. 

In this case, since the volumetric flow rate is conserved along the tree, every value of the 
intermediate flow rate qij is known from the values of qN,j, j 6 [1, 2 ]. 

The existence of solutions for the problem (<$^i) is classical. Indeed, if one of the optimization 
variables goes to zero then the energy tends to a positive infinite value. It is thus possible to come 
down to minimize £(q, •) on a compact set, which yields the existence of a solution. Moreover, £/\ 
is a convex set and £(q, •) is strictly convex, this ensures the uniqueness of the minimizer for the 
problem (<^i). 

We will now write the first order optimality conditions: we denote by the minimizer of the 
problem (^i). There exists a Lagrange multiplier A G R such that 

q 2 - 

V(i,j)eB N , --^ = A. 
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Using the volume constraint, we immediately obtain the following expression for the minimizer £*: 



The following proposition summarizes the conclusions for the first case. 
Proposition 4. The problem (&i) has a unique solution given by 



V(i,j)€B N , ^. = (A-1)=— 



(i,j)eB N W,3\ 

3.2 2nd Case: flow rate at the inlet and the pressures at the outlets are known 

The volumetric flow rate $ > which enters the overall tree through the root branch and the 
pressures at the outlets of the tree (i.e. the vector p, p £ (M?j_) 2 xl ) are assumed to be fixed. 

According to Proposition [3l the constraint q = (Mjv(^)) -1 6jy comes down to impose an in- 
compressible and Poiseuille-like flow through the tree. 

Now, let us define the set of admissible pairs (q, £) as 

= {(g,0 | £ € and q = (Mn^))' 1 b N } . (19) 
The resultant optimization problem writes 

<p.\{ min£(q,0 

We claim that generically with respect to the vector p, the problem (^2) has no solution. Never- 
theless, in the degenerate situation where all the pressures at the outlet of the tree are equal, the 
problem (^2) has a unique solution. 

We will start with the case where (^2) has a solution. The following theorem gives a charac- 
terization of this situation. 

Theorem 1. The problem (^2) has a solution if, and only if there exists tt* G K such that 

P = 7T*U N , 

where txjy € R 2JVxl denotes the unit vector (1, 1) T . Furthermore, 

( N 2 
min £ (q,£) = r <I> 2 I 1 + 



(q,flet A V A" 1 

This theorem emphasizes the fact that if two pressures at the outlets are different, then the 
problem (^2) has no solution. In this case, we are able to exhibit a minimizing sequence. It has 
to be noticed that the value of the infimum does not differ according to the case considered. 

Theorem 2. Assume that at least two pressures differ at the outlets of the tree, i.e. 

3j G [1, 2 N j I p N j ^ PN,j+i- (21) 

Then, 
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1. the problem (^2) has no solution, 



2. a minimizing sequence (£ £ ,q 6 ) £ >o is given by 

^• = E ^ ' V(»,i) GBjv\IIo-,(iV,l), ( 22 ) 

wif/i e > sufficiently small to respect the constraint: V(i,j) € Sat, • > 0, 
5. one /ias 

/ TV 2 

^. £/te sequence (q £ ) £ >o converges to the vector ($,0, . . . ,0) T as e — > 0. 

The limit case u q° = (<&, 0, . . . , 0) T " corresponds to the case where the fluid exits the tree 
only by the outlet denoted by (N,l). The use of some symmetry properties in the structure of 
the objective function £ yields that we would easily get another minimizer by choosing any other 
outlet as main exit of the fluid, which would correspond to the closure of every path linking the 
root branch to the outlet except IIo-»-(jv,i)> f° r an Y * £ [1,2 ]. Note that the symmetry property 
does not hold for the pressures at outlets. However, when we take the limit in the minimizing 
sequence, since the total flow in the tree is imposed, the pressure drop between the root and the 
outlet that remains open is the same whatever the outlet chosen. Thus, the pressure at root will 
depend on the pressure imposed at the outlet. Hence, introducing the matrix Mjy(£) that does not 
depend on the root pressure is well adapted. 

The optimal shape for the case with each pressures equal at outlets is obviously unstable rela- 
tively to the pressures: a slight change in one of them will shift the optimal shape from a tree-like 
structure to a pipe-like structure. From the modeling point of view, it means that a tree-like struc- 
ture verifying the hypothesis of this section should not be encountered in nature as a result of an 
evolution process such as natural selection. Indeed a perfect adjustment of pressure at outlets is 
very difficult to assure and maintain. 

The proofs of the theorems will be decomposed into two steps, whose details can be found in 
the next section: 

1. determination of a lower bound for the total viscous dissipated energy, 

2. construction of the minimizing sequence (q e ,£ e ) e >o- 

4 Proofs of theorems Q] and [2] (2nd case) 

4.1 Determination of a lower bound for the total viscous dissipated energy 

We will now focus on an auxiliary optimization problem whose resolution is useful in the proof of 
Theorem [2J Let q G r(2 n+1 -2)xi and £ £ ( R ^)(2 JV + 1 -2)xi be the two vectors 



q = (91,1,91,2, • • - ,9iV,l, • • • 3 Qn,2 n ) and £ = (6,1,6,2, • • • ,&V,1, 
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Notice that, in the rest of the paper and besides Section [4.11 the notation q points out the vector 
composed only of the flow rates at the outlets of the tree. 

Let ro > 0, $ > and A > 1. Let be the finite-dimensional constrained optimization 
problem 



where £(q, £) is defined by 



1 (q,€) e^i x^ 2 , v ; 



£(q,d)=r t> 2 + ^ (24) 



and the sets of constraints are 



i=i 



i G (M ^)(2«+ 1 -2)xl I £ = A — 1 



Proposition 5. T/ie "problem has a solution (</*,£*) verifying necessarily 

V(,; )eBl ,,^ = ^2j. (25) 
Hi A " 1 

Moreover, the value of the minimum is equal to r$<& 2 ^1 + ^-j- 

Note that there is no reason to have uniqueness of the solution of Problem (J2). 

Proof. The positivity of £(■, ■) grants the existence of the lower bound inf{£(q, £), (q, £) G x ^2}. 
Since the constraints are uncoupled, one has 

inf £(<?,£)= inf inf £(<?,£). (26) 
^ € ^2 being fixed, let us first focus on the optimization problem 

w{£% q,€> <27) 

Since the energy £(•,£) is obviously continuous, coercive and strictly convex, on the convex closed 
set ^x-i t ne problem has a unique solution q(£). 

By virtue of Kuhn- Tucker's theorem, there exists a real Lagrange multiplier [i\ such that 

d£ 2r <zi,i(0 , 3£ 2r <Zi, 2 (0 

a = 1 = ^ and a = 1 = Vi- 

<%,i £1,1 091,2 6,2 
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Since <7i,i(£) + <Zi,2(£) = one has 



1,1 



$ and 9l 2 (£) 



1,2 



Similarly, 



6,1 + 6,2 

Vie[l,JV], Vje[l,2i, 



6,1 + 6,2 

6,j 



<I>. 



Sfc=i 6,& 

Let us now focus on the minimization of the function £ G ^2 l— ► £(q(£)i£)- From (|28p . one gets 



(28) 



= r $ 2 ^1 + 
Let us introduce the change of variables 

2' 



6j 

(i,i)eSjv (ELi6,fc 

1 



/ 



6,1 + 6,2 



+ • • • + 



sr 2N t 

2^=l6v,i, 



Vi = i € [1,2 

i=i 



(29) 



Then, our optimization problem becomes 



& = a-i. 



mm 

N 



(30) 



^ i=l 



A direct application of Kuhn- Tucker's theorem yields that the unique minimizer of this problem 
is y* = (^p, • • • , ^r) and that the minimum is equal to ro^ 2 ( 1 + -^n) • 



From Equation ([28]) we know that q* ■ = jir<&, and we deduce equation ([25]) : jjr- = -rzT- 



Consequently, the minimum of £ in 1£i X ^2 is ro$ 2 ^1 + 
proposition follows. 



and the conclusion of this 

□ 



Remark 2. The expression of the necessarily first-order optimality conditions in the proof of Propo- 
sition^ leads easily to conclude that Problem (£}) has an infinite number of minimizer s. 

Moreover, an interesting minimizing sequence (q e , £ £ ) £ >o in ^1 x of £ for our problem is 
given by 



Ve > 0, 



V(i, j) e IIo->(Ar,i)> 9fj = $ and 3j - iv 
V(i,j) G SAr\n ^ (Aril) , = and ^ = e. 



A-l /2 JV + 1 -2 



(31) 



Physically, it corresponds to the case in which the fluid exits the tree only by the outlet denoted by 
(AM). 
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4.2 Proof of Theorem [T] 



As a preliminary to the proof, let us deal with the qualification issue for Problem (^2)- In fact, any 
elements of the set verifies the constraint qualifications. Indeed, thanks to Proposition O one 
can write q = M]\r(^)~ 1 b]\f and furthermore, it is obvious that any element of the set srf\ verifies 
the constraint qualifications. Combining the two previous remarks, the set %y inherits from this 
property . 

Let us assume that Problem (^2) has a solution (q* , £*). In the proof of Theorem [2] presented 
in Section 14.31 it is proved without any additional assumption than those of Theorem HJ that the 
sequence (q 6 ,£ £ ) £ >o defined by ([22]) is a minimizing sequence for Problem (^2), whatever values 
of pressures at the outlets. Hence it follows that 



£(q*,e) = ^£(q e ,e)=r ^ 2 1 



£->0 



N 2 



A- 1 



Now, recall that, by Proposition^ any solution of Problem (=2) realizes the minimum ro$ 2 
and verifies the necessary optimality conditions ()25|) . Since £(<?*,£*) = ro^ 2 ( 1 + 



1 + 

N 2 
A-l 



N 2 
A-l 

and 



(<f*,£*) £ C ^1 x ^2, the pair (g*,£*) belongs to the set of minimizers of Problem 
and verifies therefore, the necessary optimality conditions (|25p . namely 



V(i, j) G Sat, g. 



N<S> 



(32) 



Let us denote as previously by po the pressure at the inlet of the tree. According to Poiseuille's 
law and reasoning by induction, we obtain 



Po 
Po 



Pi,i 

Pl,2 



jVjg 
A-l 
N<5> 
A-l 



Pl,l = Pl,2 



Pi,j — Pi+l,2j-l 
Pi,j — Pi+l,2j 



In particular, it implies that 



A-l 
N<S> 
A-l 



Pi+l,2j-l =Pi+l,2j,V(i,j) 6 [J 0jv,i- 

ie[i,JV-i] 



PN,1 = PN,2 



PN,2 N -1 — PN,2 N - 



Conversely, let us prove that, if there exists tt* G R such that p = tt*un, then, Problem (^2) has 
a solution. In this case, we are able to exhibit an element (<?*,£*) which belongs to the admissible 
set %y. For instance, let us choose 



V(i,j)eB N , £*, 



A-l 



Then, it is easy to make the calculation of the flow throughout the tree. Indeed 




Q*n,i = 1*N,2 since PN,i = PN,2 = tt* and Cn ,1 = £. 



AT,2- 
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Iterating this reasoning proves that q* is proportional to un and finally, by induction, we obtain 

V(i, j) G B N , q*j = — . 

A direct calculation shows hence that 

2 A iV 2 



£(q*,r) = r <£ 2 1 + 



A-l 



4.3 Proof of Theorem [2] 



The fact that Problem (^2) has no solution when assumption (12ip holds is a direct consequence of 
Theorem [TJ In this section, we prove that the sequence exhibited in (j22j) is a minimizing sequence 
for Problem (^2), in other words, 

fim£ (</(rU £ ) = r $ 2 (l + , (33) 

where q(£ e ) is the unique solution of the linear system M^{$ i £ )q = (see Proposition [3]) . Indeed, 
since x ? 2 and since the sequence (q(£ e ), $ £ ) £ >o has been constructed so that its elements 

belong to one has for e sufficiently small, 

/ AT 2 
£(q(e),£ £ )>r $ 2 1 + 



A-l 



Hence, proving (|33p will yield at the same time that (q(£ £ ), £ e )e>o is a minimizing sequence for 
Problem (£^2) and that the infimum is equal to r^ 2 ( 1 + 



Let us denote for the sake of clarity, q(£ £ ) by q e , An(£ £ ) by An(e) and Mn(£ £ ) by Mjv(e). Now, 
recall that, according to Proposition [3l the vector q £ is completely characterized by the system 

f (q £ ,A N (e) Vi ) = -(p,Vi), Vi G [1,2^ - 1] 

\ (q e ,itAr) =*. 1 J 

Let e > 0. The matrix Aj^(e) admits the decomposition 

An(s) = r (u N + \A\ + — J +i _ 2 -^j (35) 
where Un, AL and A^ are three matrices independent of e with same size as Ajy(e), such that 



The first part, Un is the matrix of size 2^ x 2 N whose coefficients are all equal to 1. This 
term comes from the root branch that add a resistance ro to the resistance of any pathway 
in the tree. 

The second part, ^A l N , corresponds to pathways consisting only in branches such that 
= £ (i.e. whose diameter will go to with e), namely all pathways in the tree except 
those included in the pathway going from the root to the outlet (N, 1). A^ is more precisely 
characterized in Lemma Q] below. 
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• The third part — N)e n corres P on ds *° the pathways that stay open when e goes to 
0. 

Let us define A^{e) = eAn{e). Then, we have 

lim A N {e) = r A l N . (36) 

Therefore, the vector q £ is completely characterized by the system 

(tf,A N (e)vi) = -e(p,Vi), Mi e [1,3] 



(q £ ,u N ) = $ (37) 

Thus it is necessary to study more closely the matrix A]^ that will give the behaviour of q £ when 
e goes to 0. 



vJV-2 

whose branches have all the same resistance equal to 1. 



Lemma 1. A X N = (0) © @p=Q B p where B p is a resistance matrix of a tree of p + 1 generations 



Proof. Prom its definition, A^ corresponds to the resistance matrix of a tree with N + l generations 
where the resistances of the branches on the path from root to outlet (N, 1) are all and the 
resistance of the other branches are all 1. □ 



o/ \i 




2i =(o) © (i) © 



2 1 
1 2 



A 1 



( \ 

10 

2 1 

\ 1 2 / 



Figure 2: Example of the matrix A l N for N = 2. The numbers next to the branches represent their 
resistance. The matrix A l N is the direct sum of the resistance matrices of the subtrees built up by 
branches with non zero resistance. 



where 



Now, system (f3"T|) is equivalent to 

M N {e)q £ = b N {e) 
( (A N {e) Vi y \ 



M N (e) 



\ 



(^Ar(e)f 2 JV_ 1 ) T 

U N T J 



and bjv(e) 



-e(p,v 2N _ 1 ) 



(38) 
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One has 



with 



lim Mn(e) = Mn and lim bjv( e ) = ^>N 

e— >0 e— >0 



M N 



\ 



r {A l N v 2 N_ 1 ) 



and 6jv 



( o \ 


V * J 



Proof. Let us show that the family y(A^jVi) i=1 2 n -i, i^jv ) is linearly independent. We consider 
some real numbers (Aj) i=1 ...2^-1 and ^ such that 



Lemma 2. The matrix is invertible. 



2«-l 



HU N + ^2 \iA l N Vi = 



(39) 



i=l 



N-2 



Let us recall that, from Lemma[U = (0)© B p where B p is a resistance matrix and is thus 

p=0 

invertible. Consequently, the range of A l N is Span (e 2 , . . . , e 2 jv), which implies that, in equation 
([39]) the component along e\ is reduced to u = 0. 

Now, the restriction of ^4]y on its range is a one-to-one correspondence as soon as the projection 
of the family (i?i)j = i ...2^-1 on the range of is a basis of the range of Aj^, which is clearly true 
in our case. Then, the family {A l N v{) i=l ___ 2 N -i ls linearly independent and all Aj's are 0. □ 



Since the map A 1— > A 1 defined on the set of invertible matrices is continuous, the unique solu- 
tion of Misr(e)q £ = fojv converges to the unique solution of the system Mjy (lim q £ ) = by virtue 

\e— >0 / 

ofLemma[2j Furthermore, we know that MnE\ = {r^v-^ A l N e±,--- , r$v 2 N_ A l N e±,UN T e±) = 
e 2 N because A l N = A l N and Ker(Ajy) = Span(ei). Finally, 



lim q e = M Ar 1 bj V = $M N 1 e 2 N = $ei. 



£-5-0 



(40) 



Then, since we have 6jv = £7 + ^e 2 N, with 7 € Span (ei, . . . , e 2 N_ 1 ) , we have q e = M^{e) 1 bjy 
eM N (ey 1 -f + $ei and 

£(q £ ,e) = (q £ ) T A N (e)q £ 

= ^ 7 T {Mn^Y + $ei T ) A N (e) (eM^'S + $ei) 

+e$ei T A A r(e)M Ar (e)" 1 7 + ^ 2 e 1 T A N (e)e 1 



Since lim (eA^r(e)) = A\, and limMjv(e) = M^r, then one successively has 
e^o £^0 
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eW (M N {e)- l ) T A^MNie)-^ — > 



e$7 T (M J v(£)- 1 j T ^ 7V (e)ei ^ 

e$ei T A A r(e)M 7 v(e)- 1 7 — ► 



because ^4]y e i = 0, e-i T A X N = and ei T yliv(e)ei = tq + j- 
Hence, 



£->0 



/ TV 2 

lim£(qr e ,£ e ) =r ^ 2 ( + 1 ) =m ' 



which concludes the proof. 



5 Case of a fluid driven by Navier- Stokes equations: some numer- 
ical results 

As stated in the introduction, we study in this section an infinite dimensional optimization problem. 
Now the unknown is the whole shape of the dyadic tree and we do not anymore neglect the 
connections at bifurcations nor do we constrain the branches to remain cylindrical. Moreover, we 
will now use the full non linear Navier-Stokes equations and still minimize the dissipated energy 
of the fluid. Our goal is to compare the shapes obtained in this case by numerical means to that 
obtained theoretically for Poiseuille's laws in the previous sections. Nevertheless, we will not be 
able to catch the behaviour of the case with identical pressures imposed at each outlets since it is 
unstable. It is not possible in our numerical simulations to make the pressure exactly equal because 
of the rounding errors and the local meshes. Only a dedicated algorithm that would include the 
symmetry and prevent the closing of branches would be able to catch this behaviour. Consequently, 
from the theoretical study, we expect to observe the closing of every branches except one. 

We focus here on the 3D case but the 2D case can be easily adapted from this 3D study and 
numerical results will be presented both in 2D and 3D. 

The partial differential equations describing the behaviour of the fluid and the boundary con- 
ditions is now more general. In particular, this model is a convenient choice for the modeling of 
the bronchial tree (see for instance [TH1 [T9l l20l 122] ). 

Let us recall some classical definitions in incompressible Fluid Mechanics. 

Definition 4. Let u be a smooth vector field 0/M 3 (for instance C 1 ). We define 
1. the stretching tensor of u (symmetric part of the gradient tensor): 

(I ( dui du.j 



v 2 \dxj dxi j y i<jj<3 
2. the doubly contracted product of two stretching tensors e(u) and e(v) 

3 3 



1=1 J=l 
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3. the stress tensor of (u,p), where u is a smooth vector field o/R 3 representing the velocity of 
a fluid whose viscosity is fi > and p a function representing the pressure defined on R 3 : 

a(u,p) = -pl 3 + 2fie(u), 

where I3 is the identity tensor o/R 3 . 

5.1 The shape optimization problem 

We now give some precisions on the frame of our study. Let Q be a generic three dimensional 
K-shaped domain. We will denote by <9f2 the boundary of f2. In the sequel, we will assume that 
the inlet E of il is a disk and that the outlet S consists in two identical disks. We decompose the 
boundary of as the disjoint union d£l = E U V U S. T, the lateral boundary, is the main unknown 
or the shape we want to determine. In the following, we will thus consider that E and S are fixed 
part of the boundary. 

We assume that a fluid driven by the Navier-Stokes equations flows in f2. In particular, this 
model is convenient to represent the upper part of the bronchial tree since the velocity of the air 
at the beginning of the trachea is important enough to consider that the flow is inertial or even 
turbulent (the Reynolds number in trachea is about 1000 at rest). Then, the velocity u : £1 — y R 3 
and the pressure p : f2 — > R, are solutions of the Navier-Stokes system 

—fiAu + Vp + pu ■ Vu = x G Q 

v • u = x £ n 

< u = uq x e E (42) 

u = x G r 

— pn + 2[ie(u)n = h x G S, 

where 

• n denotes the outward-pointing unit normal vector at a given point of the boundary dVt, 

• uq is a parabolic velocity profile (i.e. a Poiseuille's flow is imposed at the inlet E), that is 

u = T (0,0, c(r 2 -R 2 )), 
where c is a negative constant so that the flow is ingoing, and R the radius of the inlet. 
Remark 3. Let us clarify the choice of the boundary conditions for this model. 

• The condition u = onT is the so-called no-slip boundary condition and means that the fluid 
"sticks" to the wall. 

• Practically speaking, we will impose h = —pon at the outlet S. This condition comes more or 
less from the assumption that the pressure is the sole force acting on S and that it is known 
and equal to pq. Drawing a parallel with the modeling of the bronchial tree, these conditions 
simulate muscles applying the same force per unit of surface for each outlet (—pQn ) and thus 
using the same energy for each (see JMEj). Similar boundary conditions are more detailed in 
|3 chapter 5] and in J$ 0/. 
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Notice that the classical theory of Navier-Stokes equations (see I32|) gives information on 
the existence and uniqueness of the solution of System (I42D : 



Theorem 3. Let Q, be a bounded Lipschitz domain o/R 3 . Let us assume that uq belongs to the 
Sobolev space (H 3 / 2 (E)) 3 and h belongs to (H 1 ^ 2 (S)) 3 . There exists A > such that if the viscosity 
[A is larger than A, then Problem (JW has a unique solution (u,p) G (i? 1 (il)) 3 x L 2 (Q). 

Let us now introduce the shape optimization problem we want to solve. The objective functional 
J(Q) is the energy dissipated by the fluid, i.e. 

J (SI) =2(i f |e(u)| 2 dx (43) 
Jn 

To make the statement of the optimization problem precise, we need to define a class of admissible 
shapes. As classically in shape optimization, we fix the measure Vo > of f2. Let us introduce 

Oy = {0 bounded and simply connected domain in R 3 containing E and S, such that meas(f2) = Vo} 

(44) 

The shape optimization problem writes 



min J(p,) 

n e o Vo . 



(45) 



The question of knowing if Problem (1451) has or not a solution is still an open problem. Nevertheless, 
it is possible to show that a problem, very close to Problem (J45]) but a little bit constrained, has a 
solution. 

Restricting the set of admissible shapes is a very common approach in shape optimization, 
since these problems are often ill-posed (see for instance [HD3]). A very close existence theorem is 
announced in [H] and proved in [15] , considering instead of the set Oy , a set of domains verifying 
an e-cone property, which yields some kind of uniform regularity (see [U [T3"l IS] for some reminders 
on the e-cone property and its consequences on the existence of optimal shapes) . 

5.2 Computation of the shape derivative of J 

The classical algorithm we will recall in Appendix [B] is based on the use of the well known shape 
derivative. We recall here the expression of such a derivative. The details of its calculation are 
given in [15J. 

It is possible to define an adjoint state for System (|42j) that is useful to write the shape derivative 
in a very usual form (see |13[ Theorem 5.9.2]). Let (v,q) be solution (in the case where it exists) 
of the linearized Navier-Stokes system 



' -fiAv + (Vu) T v - (Vv) u + Vq = -2[xAu x E 9, 

V • v = x £ O 

v = x £ EUF 

—qn + 2fie(v)n + (u ■ n)v — 4/j,e(u)n = x € S. 

The following existence and uniqueness result is established in [HI Proposition 3.1]. 



(46) 



Proposition 6. Let Q be a bounded Lipschitz domain ofM 3 . There exists A > such that, if the 
viscosity \x is larger than A, then the problem ( f^6| ) has a unique solution (v,q). Moreover, this 
solution belongs at least to (i^T 1 (Q)) 3 x L 2 (£2). 
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Notice that the restriction in the size of the viscosity fi ensures the existence, the uniqueness 
and a sufficient regularity of u, solution of the Navier-Stokes equation 02] (see theorem [3]) . 

We are now able to define the derivative of J with respect to the domain. Let us consider a 
regular vector field V : R 3 — > R 3 with compact support which does not meet neither E nor S. For 
small i, we define = (/ + tV)Q, the image of by a perturbation of identity and f(t) = J(Qt)- 
We recall that the shape derivative of J at with respect to V is /'(0). We will denote it by 
dJ(£l; V). To compute it, we first need to compute the derivative of the state equation. We use 
here the classical results on shape derivatives as in [13], [24], [31]. The derivative of (u,p) is the 
solution of the linear system 



-/iAu' + u' Vu + u- Vu' + Vp' = x € fl 

V • u' = ie!l 

u' = x G E 

u = — V ■ n) 



(47) 



dn 

p'n + 2fie(u')n = 



xer 



where u is, under assumption that fi is large enough, the unique solution of (|42p . Similar formula in 
the context of Shape Optimization in Fluid Mechanics are established for instance in 
Now, we have (see [13], [51] ) 



(48) 



dJ(ft;V) = 4/i / e(«) : s(u')dx + 2/x / \e(u)\ 2 (V ■ n)d& 
Jn Jr 



Using the adjoint state (|46|) . it is possible to rewrite dJ(f2, V) in a more workable form (see |15t 
Proposition 3.2]): 

d,J(n; V) = 2fx J (e(u) : e{v) - \e{u)\ 2 ) (V ■ n)ds. (49) 

We will now introduce the shape gradient of the criterion J as a functional defined on the boundary 
such that d J(ft; V) = j r V J(Q) ■ Vds, in other words 



VJ(fi) = 2/i (e(u) : e(v) - \e(u)\ 2 ) n. 



(50) 



Thanks to this expression of the shape gradient, we are now able to implement an augmented 
Lagrangian algorithm in order to solve Problem (|45h . The algorithm is classical and is described 
in Appendix [Bj 



5.3 Numerical results 

In this section, we present 2D and 3D numerical simulations using the augmented Lagrangian 
algorithm described in Appendix[B]based on the shape gradient (|50p . They have been implemented 
with the script of the software Comsol Multiphysics. The 2D algorithm is an easy adaptation of 
the 3D algorithm. 

The Navier-Stokes systems and the adjoint state are solved with a direct finite elements method 
using Lagrange elements PI for pressures and P2 for velocities. The displacements of the mesh 
are PI. The multifrontal package (UMFPACK) is used to solve the resultant linear system and 
a modified Newton like method is used to treat the non linear term. At each iteration k, each 
node of the geometry is perturbed by the discretized operator (/ + £fcdfc)(J7fc) (classical Arbitrary 
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Lagrangian Eulerian method) and it is necessary to remesh the geometry when the displacement 
of the mesh becomes too large. 

The Reynolds number in the following computations is 100 which is not only large enough to 
observe inertial effects near the bifurcation, but small enough to let the computations run in a 
reasonable amount of time. The effect of inertia on flow distribution in the bifurcation can be seen 
on Figure [3] where the image on the left represents a non inertial flow (Reynolds 0) and the image 
on the right an inertial flow with Reynolds 100. 



Figure 3: Velocity norm (jet colors mapping) in a 2D bifurcation for a non inertial flow with 
Reynolds (left) and for an inertial flow with Reynolds 100 (right). Inertia alters the flow mostly 
around the bifurcation. 

We normalized the applications J(-) and V(-) in such a way that J(f2o) = 1 an d V(£Iq) = 1. 

5.3.1 1st Case: Neumann conditions are imposed at the inlet and Dirichlet conditions 
are imposed at the outlets 

In this section, the initial geometry 0,q is a bifurcation whose branches are identical except for the 
length of the mother branch which is 1.5 the length of the daughter branches, see Figure 0] (2D, left 
image) and Figured] (3D, left image). The length of the mother branch has been chosen longer to 
avoid that the flow near the bifurcation interacts too much with the inlet flow since the algorithm 
tends to shorten the mother branch. Practically, a mother branch too short leads to convergence 
problem. The angle between two nearby branches is tt/3. 

The boundary conditions correspond to that of Case 13 . ll adapted to the non- linear regime. Thus 
a parabolic velocity profile is imposed at each outlets (Dirichlet condition) and Neumann boundary 
conditions a(p, u)n = —pon are imposed at the inlet, which is "formally" interpreted as a pressure 
constraint since viscous effects are negligible. The pressure imposed at the inlet is 0. 

Figures [5] (2D) and [7] (3D) show the convergence of the different quantities of the problem: the 
Lagrange multiplier (upper left), the Lagrangian function (upper right), the viscous energy (lower 
left) and the volume (lower right). Except the Lagrange multiplier, these quantities have been 
normalized. As expected, the curves are oscillating around their asymptotic values since verifying 
the volume constraint and minimizing the criterion are advantaged one after the other by the 
algorithm. This "advantage" is controlled by the Lagrange multiplier. Hence when it is too large 
the volume constraint is stronger and the algorithm authorized an increase of the criterion; if it is 
too small then the algorithm let loose the constraint and decrease the criterion. The amplitude of 
these oscillations decreases with the iterations. 

In 2D, the mesh consists in 9296 triangle elements and no remeshing was needed to ensure 
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convergence. Convergence was reached for 1000 iterations of the algorithm in about 6 hours using 
2 cores of a Xeon processor (2.33 GHz). The viscous energy dissipated in the bifurcation, which 
is the quantity minimized by our algorithm, has been reduced by 20.4% in the final geometry in 
comparison with its value in the initial geometry. Volume precision at convergence was smaller 
than 2 x 10~ 6 . 

In 3D, the mesh consists in 17959 tetrahedral elements and no remeshing was needed either. 
Convergence was reached for 1700 steps. Computation time was about 46 hours using 2 cores of a 
Xeon processor (2.33 GHz). The viscous energy dissipated in the bifurcation has been reduced by 
22.5% between the initial and final geometry. The volume constraint at the end of the algorithm 
was smaller than 3 x 10~ 5 . 




Figure 4: 2D Case. Left: initial geometry Qq, the inlet is the blue branch while the outlets are the 
red branches. Right: final shape reached by the algorithm, the colors represent the norm of the 
velocity u (increasing with the warmth of colors). 



5.3.2 2nd case: Dirichlet conditions are imposed at the inlet and Neumann conditions 
are imposed at the outlets 

In this section, the initial geometry f^o is a bifurcation whose branches are all identical, see Figure 
[8] (2D, left image) and [11] (3D, left image). The angle between two nearby branches is n/3. 

The boundary conditions correspond to that of C ase 13 . 2l adapted to the non- linear regime. Thus 
a parabolic velocity profile is imposed at the inlet (Dirichlet condition) and Neumann boundary 
conditions a(p, u)n = —pon are imposed at the outlets. The pressures imposed at the two outlets 
are slightly different: Pa on one outlet and 2 x 10 -4 Pa on the other. 

As for the first case, the curves plotted on figure [TUl (2D) and[12](3D) represent the convergence 
of the different quantities of the problem: the Lagrange multiplier (upper left), the Lagrangian 
function (upper right), the viscous energy (lower left) and the volume (lower right). Except the 
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Figure 5: 2D case. Convergence curves of the Lagrange multiplier the Lagrangian £5(0^.,^), 
the viscous energy J(fifc) and the volume V(f2fc). The x-coordinates represent the iterations number. 




Figure 6: 3D case. Left: initial geometry Qq, the inlet is the blue branch while the outlets are the 
red branches. Right: final shape reached by the algorithm. 



Lagrange multiplier, these quantities are still normalized. As before, they are oscillating around 
their convergence value (see 1st case). 

In 2D, the mesh consists in 8624 triangle elements and no remeshing was needed. The conver- 
gence was achieved in 2800 steps for a computing time of 8 hours on 2 cores of a Xeon processor 
(2.33 GHz). In the final geometry the viscous dissipated energy is reduced by 54.5% relatively to 
the initial bifurcation. The precision of the volume constraint is smaller than 7 x 10 -5 . 

In 3D, convergence was achieved in 7000 steps. A remeshing was performed at the step 4500 
(represented by the dotted vertical line on figure fTTj) . The initial mesh consisted in 10044 tetrahedral 
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Figure 7: 3D case. Convergence curves of the Lagrange multiplier the Lagrangian £^(0^, /ij.), 
the viscous energy J(Qk) and the volume V(f2fc). The x-coordinates represent the iterations number. 



elements and the second mesh was finer with 37167 tetrahedral elements. The computation was 
much slower after the remeshing because of the increased number of elements. The total time 
needed by the algorithm to converge was of 259 hours on 2 cores of a Xeon processor (2.33 GHz). 
The viscous energy dissipated in the final geometry is reduced by 53.6% relatively to the viscous 
energy dissipated in the initial geometry. The volume constraint precision at convergence was 
smaller than 7 x 10 -5 . 

Similarly than the result found theoretically at low regime in Section 13. 2\ we observe the closing 
of one branch. This indicates that this result should be more generally true and that it can probably 
be extended to non linear regime. A zoom of the 2D bifurcation is plotted on Figure El the arrows 
represents the normalized velocities and a region of fluid recirculation appears at the beginning of 
the closing branch. The colors on Figure [9] represents the pressure which is almost constant along 
the branch that is closing (up). This induces a very small flow inside it (the flow and the branch 
diameter are decreasing together when the optimization algorithm is progressing). 

6 Conclusion 

In this work, we show that fluid boundary conditions play an important role for the determination 
of the optimal shape of a tree in term of viscous dissipation. Indeed, the optimal shape associated to 
pressure conditions at leaves is a simple pipe while the optimal shape associated to flow conditions 
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Figure 8: 2D case. Left: initial geometry Qq, the inlet is the blue branch while the outlets are the 
red branches. Right: final shape reached by the algorithm, the color represents the norm of the 
velocity u (increasing with the warmth of colors). 




Figure 9: Zoom on the bifurcation, the color represents the pressure and the arrows the directions 
of the velocities (normalized vector field). A recirculation zone is present in the bifurcation and 
the pressure is almost constant in the closing branch (up) inducing a very small flow in it. 



at leaves is a tree. Moreover we have shown that these results hold for both Poiseuille's regime 
and a regime with a Reynolds number 100, which, although moderate, is large enough to exhibit 
inertial effects near the bifurcation. 

Boundary conditions can be seen as constraints in the system and a pressure constraint is very 
different of a flow constraint. Let us consider a pipe in which flows a fluid in Poiseuille's regime. If 
we call R the hydrodynamic resistance of the pipe, the dissipated energy is J = RQ 2 = (Ap) 2 /R 
(<I> is the flow, Ap is the pressure drop). Then minimizing J relatively to the geometry of the pipe 
depends on the constraint: 

Situation (i): if the flow $ is given then minimizing J is equivalent to minimizing the resistance 
R and the minimizer corresponds to Ap = (pipe radius goes to infinity). 

Situation (ii): if the pressure drop Ap is given then minimizing J is equivalent to maximizing the 
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Figure 10: 2D case. Convergence curves of the Lagrange multiplier /i^, the Lagrangian /^(O^, 
the viscous energy J(fJfc) and the volume V(f2fc). The x-coordinates represent the iterations number. 




Figure 11: Left: initial geometry £1q, the inlet is the blue branch while the outlets are the red 
branches. Right: final shape reached by the algorithm. 



resistance R and the minimizer corresponds to <3? = (no flow in the pipe, the pipe radius 
goes to 0). 

Simply speaking, the case of the tree we studied in this paper is an extrapolation of these two 
points depending on the conditions imposed at exits. When flows are imposed at outlets then each 
outlet is in Situation (i) and the optimal geometry is the one that tries to open the branches (in 
the limit of the constraint on the volume). When pressures are imposed, then all outlets are in 
Situation (ii) except one which is in Situation (i) since we imposed a non zero flow in the root of 
the tree and this flow, by conservation, has to get out of the tree by at least one branch. Moreover, 
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Figure 12: Convergence curves of the Lagrange multiplier fi^, the Lagrangian C^ilk-, ^k)-, the 
viscous energy J(Ofc) and the volume V(Qk)- The x-coordinates represent the iterations number. 
The dashed line represents the step 4500 at which a remeshing was done. 



in this last situation, our results show that it is "better" to close a branch and to use its volume to 
widen another branch. This implies that the other branch can distribute a larger amount of flow 
while dissipating less energy than two separate branches. Indeed, the viscous effects are larger near 
the walls and one wide branch has less walls than two smaller branches. When flow is imposed in a 
branch, this phenomena is compensated by the fact that reducing the branch radius increases the 
viscous effects by increasing the velocity gradients in the flow. Consequently, the optimal geometry 
is a compromise. These phenomena are probably true whenever the inertial effects are present or 
not, as our numerical simulations confirmed partially (Reynolds 100). 

Finally, in term of modeling, our work shows that very different optimal structures can be 
obtained by boundary conditions adjustment (tree or pipe). For organs, we could imagine that 
depending on the organ function, optimization has been made through evolution by minimizing 
identical costs with different boundary conditions. 
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A Model proofs 



A.l Proof of Proposition [T] 

The linearity of the relation between p and q comes from Poiseuille's law. Hence it is sufficient to 
compute the pressure vectors p related to the elements of canonical basis of M 2 x 1 . 

Let us begin with q = (1,0,... ,0) T G R 2 xl . It corresponds to the case in which the fluid 
exits the tree only by the outlet denoted by (N, 1). In other terms, the fluid flows through the tree 
using the path IIo^^,!)- By conservation, the volumetric flow rate which enters the tree through 
the root branch is exactly $ = 1, so that the pressure p\ at the outlet of the root branch is Pq — vq. 
As there is no flow in the right-hand subtree stemming from the root branch, the pressure at its 
outlets (whose couples of indexes are between (N, 2 N ~ l + 1) and (N,2 N )) is the same as at the 
outlet of the root branch, namely po — r o- Similarly, the pressure p\ \ at the outlet of the branch 

denoted by (1, 1) is po — (ro + r% t i) = po — tq (l + • Following this approach recursively, one 

finds pressures at the outlet of the branches of the path IIo^^^) denoted by (i, 1) with i G [1, NJ 

to be po - (r + EfrQeno-wxjW r k,l) . whi ch writes again p - r (l + £(fc,0ello- W i)«) 5j) • 
Therefore, the pressure pnj (with j G Jl, 2 ]) at the outlet of the tree is 

PN,j = Po - r 1 + ^2 

\ (*,0eno_^(Ar l x)(JV-^,j_x) 

Applying the same reasoning to any vector q = (0, . . . , 0, 1, 0, . . . , 0) T G R 2JVxl (with 1 at the z-th 
position) yields 



/ 



P 



( 



PO-ro(l+ 7" PO-^O 1+ 



y \ (fc,Oen _ KJV , i) (iv-iy i ,o) ^ k ' 1 



1 



y (fc,0en o ^(Ar,o(^-«' i ,aW_ 1 ) ^ M 



Consequently, we obtain 

2^ / 

ViG [1,2^], PN,i=Po-r J2 < lJ 1+ Yl 

i= l \ {k,l)&i ^ {N>i) {N-n-\,j-i) 

A. 2 Proof of Proposition [3] 

The set of vectors {vi) i& ^2 N -i\ forms a basis of { , ujv}" L > the orthogonal complement of u^. Let 
i G [1, 2^ — 1]. Multiplying (fTUI) by Vi in the sense of the inner product yields 

V* G [1,2" - I], -(p, Vi ) = (A N (Z)q,Vi). (51) 

Furthermore, by conservation of the flow rate, one has (q,ttjv) = *F 

Thus, since An{£) is real symmetric, the vector q verifies the following system 

(q, A N (£)vi) = -{p, Vi ), Vi G [1, 2 N - 1] 
(q,u N } = $. 
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That explains the structure of the matrix Mjv(£). Let us now to prove that Mjv(£) is invertible. 
It amounts to prove that the set of vectors (An(£)v\, . . . , An(£)v 2 n_ i, M/v) is linearly independent. 
The fact that the set (An(£)v-l, . . . , An(£)v 2 n _ x ) is linearly independent is almost trivial, since 
An(£) is invertible by Proposition [2] and since the vectors (i^gin 2 N -1J are linearly independent. 

Now, assume that 

2^-1 

Xun + HiAN{£)Vi = with (A, fix, . . . , (JyN-x) £ ^ 2 • (53) 
i=i 

Multiplying the equation (|53p by the vector An^^^-un in the sense of the inner product yields 

2 N -1 

\{u N ,A N (£)- 1 u N )+ fiiiAN^v^ANi^UN) =0. 

1=1 

Since v { G {mat} 1 " for all i G [1, 2^ - lj, it follows that 

Vi G [1,2* - 1], (^(0«i 5 ^(0 _1 «Jv) = («<»«jv) = 0. 

Moreover, (tijy, -Ajv (£)~ 1m at) > 0. Indeed, we know that An{£) is symmetric positive definite, and 
therefore A^^) -1 inherits from these properties. Consequently, A = 0. 
Now, (1531) rewrites, 

2 N -1 

^2 ViAN(€)vi = 0. 
1=1 

By the linear independence of the set (A/vCO^iXelTi, 2-^-1]) we conclude that V£ G [1, 2^ — 1], ^ = 0. 
It proves that the vectors (An(£)v-l, . . . , An(£)v 2 n_ 1 , «jv) are linearly independent. 

B An augmented Lagrangian algorithm 

We present in this section the augmented Lagrangian algorithm used in this work to optimize the 
energy dissipated by a fluid in a dyadic tree with respect to the shape. 

One of the main interest of the augmented Lagrangian if compared to the classical Lagrangian 
is the regularization operation which generally improves the condition number of the dual function. 
Thus, a faster convergence of the maximizing sequence of the Lagrange multipliers is expected. 
However, the choice of a good parameter of augmentation can be very difficult. We will discuss 
this aspect concerning our simulations at the end of this section. 

The descent direction in the main step of this algorithm will be computed thanks to a gradient 
method, which implies to calculate at each iteration the derivative of our criterion with respect to 
the domain. 

The augmented Lagrangian associated to Problem (j4"5j) is 

£ b (n, £) = J(Q) + £G(Q) + h - (G(n)) 2 , (54) 

where £ G R is the Lagrange multiplier, associated with the volume constraint G(Cl), that is 

G{n) = meas(Sl) - V . 
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It is then easy to determine the shape derivative of £5. One has 

dC b (n-,V) = J [2fi(e(u) : e(v) - \e(u)\ 2 ) +£ + 6(meas(fi) - V )] (V.n)ds. (55) 

We give now some precisions on the second step of the augmented Lagrangian algorithm, in 
particular on the choice of the descent method. Let £l k be the domain obtained at iteration k — 1, 
its lateral boundary and \i k the associated Lagrange multiplier. Qk+l is searched as a perturbation 
of the identity. That is why we write £lk+i = (I + e kd k )(£lk), where d k is a vector field representing 
the perturbation of the mesh and a variable step. 

Since we want to implement a gradient method, a first approach would be to consider d k such 
that 

d fc | rfc =-VA(«fc,4), VfceN. (56) 

This question has been much studied (see for instance pjj HI [TUl E21 [231 [26]). I n particular, in [TU] . 
the authors study similar methods applied to image segmentation, and exhibit some situations in 
which the choice of d k as in (|56p is the worst solution from a numerical point of view. 
In this work, we chose d k such that 



which corresponds to d k solution of the equation 



d k ), 



-Ad k + d k = 
d k = 

-q — = — v£b(rifc,< 



x G £1 k 
x^EUS 

x G T k - 



(57) 



We are now able to write the algorithm of resolution of the Problem (jU 



Augmented Lagrangian algorithm for the resolution of Problem ( B5j) 

1. Initialization. Choose SIq £ E and £q G R. 
Let also fix r > and e s top- 

2. Iteration k. £k is known. 

(a) Resolution of the Navier-Stokes problem (and storage of its solution u k ) 

—fiAu k + Vpfc + u k ■ Vitfc = i€S!(. 
V • u k = x G 0^ 

u k = Uq x G -E 

it fc = 

-pfcn + fie(u k )n = —po-n x G 5. 

(b) Resolution of the adjoint state (and storage of its solution v k ) 

-fiAv k + (X7u k ) T v k - (V«fc) itfc + = -2fiAu k 
V • ^fc = 

^ = x G £ U T k 

-qkn + /J,s(vk)n + (life • n)ufc - 4/j,e(u k )n = x G 5. 
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(c) Calculation of the scalar 

/3 k = V£ b (Sl k ,n k ) ■ n 

= 2/i (e(tt fe ) : e{v k ) - \e(u k )\ 2 ) + 4 + 6 (meas (O fc ) - V ) . 

(d) Determination of the displacement of the mesh d k as the solution of the elliptic equation 

— AtZfc + dk = x G 0^ 
<9d fc 



dn 



-f3 k n x £ T k . 



(e) Determination of an e k that decreases the augmented Lagrangian. 

(f) Determination of the domain ^l k +i- Qk+i = (J + £fcdfc)(f2fc). 

(g) Reinitialization of the Lagrange multiplier: 4+1 = 4 + t (meas (fife+i) — Vq). 



3. Stopping criterion. The algorithm stops if |4+i — 4| < £stop and loops back to step [2] if the 
inequality is false. 

As pointed out upwards, the choice of the parameter 6 in the augmented Lagrangian algorithm 
can be difficult. Practically speaking, b has to be chosen neither too big nor too small. Indeed, the 
biggest is the parameter b, the best is the conditioning of the dual functional giving the constraints 
and then the convergence of the sequence of Lagrange multipliers. Nevertheless, if b is chosen too 
big, the conditioning of the primal problem min{£;,(f2, 4)> ^ € E} deteriorates and it becomes 
more difficult to solve. Thus a compromise has to be done. Practically, a lot of preliminary tests 
have to be done to find b before the algorithm could be run properly. 
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